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Abstract 

We study the spike statistics of neurons in a network with dynamically 
balanced excitation and inhibition. Our model, intended to represent a 
generic cortical column, comprises randomly connected excitatory and in- 
hibitory leaky integrate-and-fire neurons, driven by excitatory input from 
an external population. The high connectivity permits a mean-field de- 
scription in which synaptic currents can be treated as Gaussian noise, the 
mean and autocorrelation function of which are calculated self-consistently 
from the firing statistics of single model neurons. Within this descrip- 
tion, we find that the irregularity of spike trains is controlled mainly by 
the strength of the synapses relative to the difference between the firing 
threshold and the post-firing reset level of the membrane potential. For 
moderately strong synapses we find spike statistics very similar to those 
observed in primary visual cortex. 

1 Introduction 

The observed irregularity and relatively low rates of the firing of neocortical neu- 
rons suggest strongly that excitatory and inhibitory input are nearly balanced. 
Such a balance, in turn, finds an attractive explanation in the mean-field de- 
scriptions of Amit and Brunei ^ |2 El and Van Vreeswijk and Sompolinsky 
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0JOj|. In their theories, the balance does not have to be put in "by hand"; 
rather, it emerges self-consistently from the network dynamics. This success 
encourages us to study firing correlations and irregularity in models like theirs 
in greater detail. In particular, we would like to quantify the irregularity and 
identify the parameters of the network that control it. This is important be- 
cause one can not extract the signal in neuronal spike trains correctly without a 
good characterization of the noise. Indeed, an incorrect noise model can lead to 
spurious conclusions about the nature of the signal, as demonstrated by Oram 
et al ®. 

Response variability has been studied for a long time in primary visual cor- 
tex [!Ji rT77irT7iran^n^rani)lfT7irT^ini71 and elsewhere gDI El CHI ED Most, 
though not all, of these studies found rather strong irregularity. As an example, 
we consider the findings of Gershon et al |17| . In their experiments, monkeys 
were presented with flashed, stationary visual patterns for several hundred ms. 
Repeated presentations of a given stimulus evoked varying numbers of spikes in 
different trials, though the mean number (as well as the PSTH) varied system- 
atically from stimulus to stimulus. The statistical objects of interest to us here 
are the distributions of single-trial spike counts, for given fixed stimuli. Often 
one compares the data with a Poisson model of the spike trains, for which the 
count distribution P(n) — mV m /n!. This distribution has the property that 
its mean (n) = m is equal to its variance (8n 2 ) — ((n — (n)) 2 ). However, the 
experimental finding was that the measured distributions were quite generally 
wider than this: (Sn 2 ) > m. Furthermore, collecting data for many stimuli, the 
variance of the spike count was fit well by a power law function of the mean 
count: (Sn 2 ) oc m v , with y typically in the range 1.2 — 1.4, broadly consistent 
with the results of many of the other studies cited above. 

Some of this observed variance could have a simple explanation: The con- 
dition of the animal might have changed between trials, so the intrinsic rate at 
which the neuron fires might differ from trial to trial, as suggested by Tolhurst 
et al ^T]. But it is far from clear whether all the variance can be accounted 
for in this way. Moreover, there is no special reason to take a Poisson process 
as the null hypothesis, so we don't even really know how much variance we are 
trying to explain. 

In this paper, we try to address the question of how much variability, or 
more generally, what firing correlations can be expected as consequence of the 
intrinsic dynamics of cortical neuronal networks. The theories of Amit and 
Brunei and of van Vreeswijk and Sompolinsky do not permit a consistent study 
of firing correlations. The Amit-Brunel treatment assumes that the input to 
neurons is uncorrelated in time (white noise). Thus, although one can calculate 
the variability of the firing j^j, it is not self-consistent. Van Vreeswijk and 
Sompolinsky use a binary-neuron model with stochastic dynamics which makes 
it difficult, if not impossible, to study temporal correlations that might occur in 
networks of spiking neurons. Therefore, in this paper we do a complete mean- 
field theory for a network of leaky integrate-and-fire neurons, including, as self- 
consistently-determined order parameters, both firing rates and autocorrelation 
functions. A general formalism for doing this was introduced by Fulvi Mari jjj 
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Figure 1: Structure of the model network. 

and used for an all-excitatory network; here we employ it for a network with 
both excitatory and inhibitory neurons. A preliminary study of this approach 
for an all- inhibitory network was presented previously \22\ . 

2 Model and Methods 

The model network, indicated schematically in Fig.^ consists of N± excitatory 
neurons and N2 inhibitory ones. In this work we use leaky integrate-and-firc 
neurons, though the methods could be carried over directly to networks of other 
kinds of model neurons, such as conductance-based ones. They are randomly 
interconnected by synapses, both within and between populations, with the 
mean number of connections from population b to population a equal to Kb, 
independent of a. In specific calculations, we have used K\ from 400 to 6400, 
and we take K2 = K-y/A, The population sizes N a do not enter directly in the 
mean field theory, only their ratios (the connection probabilities) K a /N a . We 
have used K a /N a = 0.1 for both excitatory and inhibitory connections, implying 
Nt = AN 2 . 

We scale the synaptic strengths in the way van Vreeswijk and Sompolinsky 
did0J|S], with each nonzero synapse from population b to population a having 
the value J a b/VK~b~- The parameters J a b are taken to be of order 1, so the net 
input current to a neuron from the Kb neurons in population b connected to it 
is of order y/Kb. With this scaling, the fluctuations in this current are of order 
1. 

Similarly, we assume that the external input to any neuron is the sum of 
Kq 1 contributions from individual neurons (in the LGN, if we are thinking 
about modeling VI), each of order 1/^/Kq, so the net input is of order \[Kq. In 
our calculations, we have used Kq = K\. 
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We point out that this scaling is just for convenience in thinking about the 
problem. In the balanced asynchronous firing state, the large excitatory and 
inhibitory input currents nearly cancel, leaving a net input current of order 1. 
Thus, for this choice, both the net mean current and its typical fluctuations 
are of order 1, which is convenient for analysis. The physiologically relevant 
assumptions are only that excitatory and inhibitory inputs are separately much 
larger than their sum and that the latter is of the same order as its fluctuations. 

Our synapses are not modeled as conductances. Our synaptic strength sim- 
ply defines the amplitude of the postsynaptic current pulse produced by a single 
presynaptic spike. 

The model is formally specified by the sub-threshold equations of motion for 
the membrane potentials uf (a = 1,2, i = 1, . . . N a ): 

6=0 j=l 

together with the condition that when uf reaches the threshold 9 a , the neuron 
spikes and the membrane potential is reset to a value uf. The indices a or b = 1 
or 2 label populations: 6 = refers to the (excitatory) population providing 
the external input, a = 1 refers to the excitatory population and a = 2 to 
the inhibitory population. In ifTJl. r is the membrane time constant (taken the 
same for all neurons, for convenience), the strength of the synapse from neuron 
j in population b to neuron i in population a is denoted by J?- , and Sf(t) = 
Y] s S(t — tj b ) is the spike train of neuron j in population b. We have ignored 
transmission delays, and we take the thresholds 9 a = 1 and the reset levels uf 
equal to the rest value of the membrane potential, 0. In our calculations, the 
thresholds are given a Gaussian distribution with a standard deviation equal to 
10% of the mean. Analogous variability in other single-cell parameters (such as 
membrane time constants) could also be included in the model, but for simplicity 
we do not do so here. 

We assume that the neurons in the external input population (b = 0) fire as 
independent Poisson processes. However, the neurons in the network (b = 1,2) 
are not in general Poissonian; it is their correlations that we want to find in this 
investigation. 

Mean Field Theory: Stationary States 

We describe the mean field theory and its computational implementation first 
for the case of stationary rates. When the connectivity is large and random, as 
we will assume here, each of the three terms in the sum on b on the right-hand 
side of Q can be treated as a Gaussian random function with time-independent 
mean. The simplest case is b = 0, the external input. For simplicity, we assume 
that all No neurons in the external population fire at the same rate, ro- But 
because of the random connectivity, the net time-averaged input current they 
provide to a neuron in cortical population a can vary from neuron to neuron. 
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Assuming large, dilute connectivity (Kq 3> 1 and Kq <C No), the central limit 
theorem then implies 



{If (*)) = E 4°^ = £( J ^ + = J«oro(y/K~o + xf), (2) 

3 3 

where xf° is a Gaussian-distributed random number of unit variance. By (• • •) 
we mean a time average or, equivalently, an average over "trials" (independent 
repetitions of the Poisson processes defining the input population neurons) . We 
will generally use a bar over a quantity to indicate an average over the neuronal 
population or over the distribution of the J°j> . (Note that these two kinds of 
averages are very different things.) 

Writing the spike train Sj (t) for neuron j in the input population as 

S°(t) = r Q + SS°(t), (3) 
with (SSj(t)) = 0, we can write the fluctuations around (If ) as 

Slf(t) = E Ji°*S°(t) = J a0 g°(t) (4) 

3 

where £f (£) is white noise of power tq: 

(Zfm°(t')) = r 6(t-t') (5) 

Thus, quite generally, the input has a large mean value, of order \/Kq, plus 
Gaussian fluctuations of order 1. The fluctuations are of two kinds. One is 
constant for a given neuron, independent of time and trial and arises from the 
fact that the connectivity is random and the neurons in the input population 
have a distribution of rates. The other fluctuation is a dynamical one, with 
correlations (independent of i) reflecting the Poisson dynamics of the input 
population neurons. 

The recurrent input terms Ii b {t) also have large means and fluctuations, 
static and dynamic, of order 1, but certain features of their statistics are slightly 
different, as a systematic formal derivation [71|H] proves. Here we do not give 
the derivation, but just describe the result, which is that If b (t) can be written 



If{t) = J ab [VK b r b + B b xf + ^1-K b /N b g b (t)}, (6) 

with 
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the average rate in population b, x°: a unit-variance Gaussian random number, 
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and 

{itm\t')) = c b (t-t'). (9) 

Here Ct(t — t') is the average autocorrelation function of the firing of neurons 
in population b, 

C b (t-t') = l-'£(SS b j (t)6S b j (t')), (10) 

3 

with 

8S){t) = S)(t)-r). (11) 

Again, if is time- and trial-independent, while the noise (t) varies both in 
time within a trial and randomly from trial to trial. Note that for this model a 
correct and complete mean field theory has to include rate fluctuations, through 
(r^) 2 , and the firing correlations, given by Cb(t — t'), as well as the mean rates. 

The means of the recurrent input currents if 3 are completely analogous to 
the mean term in If , but the effective noise is different in three ways 

1. The amplitude Bf, of the static noise component (the second term) contains 

a factor of the rms rate *J (rj) 2 , not rt, as in (J2J. The same would be true 
for the static input noise (b = 0) if we allowed a distribution of rates in the 
input population. So this difference is not an essential one. It occurs only 
because we made a simplifying assumption about the input population. 
However, we are not allowed to assume that about the neurons in the 
cortical network, which will always have a distribution of rates because of 
the random connectivity. 

2. The neurons providing the source of these currents are not generally Pois- 
sonian, so their correlations appear in the statistics of the noise term. 

3. The noise terms, both static and dynamic, have a factor y/l — K^/Nb in 
front of them. This can be understood in the following way: It is the 
randomness in the synaptic connections in the network that generates 
these noise terms in the effective single-neuron problem; in general, they 
are proportional to (SJ-j) 2 , which is equal to J^ h (l — Kb/N^/Nj, in our 
model. In the limit of full connectivity, Kb = Nb, all are equal and 
there is no randomness. Therefore there is no noise, as guaranteed here 
by this factor. 

The self-consistency equations of mean field theory are simply the conditions 
that the average output statistics of the neurons, r a , (r") 2 and C a (t — t') are 
the same as those used to generate the inputs for single neurons using integrate- 
and-fire neurons with synaptic input currents given by J2J, Q and ijfJJl. 

In an equivalent formulation, the second term in © can be omitted if the 
noise terms £f b (i) have correlations equal to the unsubtracted correlation func- 
tion 

^ ot (i-iO = ^£<s)(*)s)(*')> (12) 
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instead of fTDjL For -> oo, Cl ot {t-t') -» (r?), so £? 6 (i) acquires a random 

static component of mean square value (r^) 2 . 

In still another way to do it, one can use the square of the average rate, r\ 
in place of (r b ) 2 in Eq.© for Bf, and employ noise with correlation function 

C b (t ~t') = ] r £<(S)(t) - n){S b 3 {t') - r b )). (13) 

j 

For |t-t'| -> oo, 

C b (t-t r )^(r b J -r b y='{S^. (14) 

There are now two static random parts of If h (t) : one from the B b term and 
one from the static component of the noise. Their sum is a Gaussian random 
number with standard deviation equal to B^ as given in JSJ. Thus these three 
ways of generating the input currents are all equivalent. 

The balance condition 

In a stationary, low-rate state, the mean membrane potential described by Q 
has to be stationary. If excitation dominates, we have duf/dt oc ^/Kq, implying 
a firing rate of order V ' Kq (or one limited only by the refractory period of the 
neuron). If inhibition dominates, the neuron will never fire. The only way to 
have a stationary state at a low rate (less than one spike per membrane time 
constant) is to have the excitation and inhibition nearly cancel. Then the mean 
membrane potential can lie a little below threshold, and the neuron can fire 
occasionally due to the input current fluctuations. Thus, using © and ®, we 
have 

2 

Jab\fK~bn = 0{1) (15) 

6=0 

or, up to corrections of 0(l/v^Ko), 

2 

JabTb = (16) 

6=0 

with J a b — J a b\/ Kb/ 1 Kq. These are two linear equations in the two unknowns 
r a , a — 1,2, with the solution 

2 ^ 

r a = >^[J }abJbor , (17) 

6=1 

where J -1 is the inverse of the 2x2 matrix with elements J a b, a,b = 1,2. If 
there is a stationary balanced state, the average rates of the excitatory and 
inhibitory populations are given by 117(1 (in the large- N limit). 
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This argument depends only on the rates, not on the correlations, and is 
exactly the same as that given by Amit and Brunei and by Sompolinsky and 
van Vreeswijk. 

This calculation does not say whether this state is stable, however. To deter- 
mine this, one can expand around this solution and examine the linear stability 
of the fluctuations, as done for their model by van Vreeswijk and Sompolinsky 
Here, we do not do this analytically, but rather check the stability of our 
states numerically within our algorithm. 

Numerical procedure 

For integrate-and-fire neurons in a stationary state, the mean field theory can be 
carried out analytically if a white-noise (Poisson firing) approximation is made 
013 OH- But if firing correlations are to be taken into account, it is necessary 
to resort to numerical methods. Thus we simulate single neurons driven by 
Gaussian synaptic currents, collect their firing statistics to compute the rates 
r a , rate fluctuations (5r°-) 2 and correlations C a (t — t'), and then use these to 
generate improved input current statistics. The cycle is repeated until the input 
and output statistics are consistent. This algorithm was first used by Eisfcllcr 
and Opper 23 to calculate the remanent magnetization of a mean field model 
for spin glasses. 

Explicitly, we proceed as follows. We simulate single excitatory and in- 
hibitory neurons over "trials" 100 integration timesteps long. (We will call each 
timestep a "millisecond" . We have explored using smaller timesteps and verified 
that there are no qualitative changes in the results.) We start from estimates 
of the rates given by the balance condition, which makes the net mean input 
current vanish. Then the sum of the 0(y/Kb) terms in (J2J and JJJJ vanishes, 
leaving only the rate fluctuation and noise terms. We then run 10000 trials of 
single excitatory and inhibitory neurons, selecting on each trial random values 
of i-' and £,f b (t). Since at this point we do not have any estimates of either the 
rate fluctuations (Sr^) 2 or the correlations C b {t — t'), we use r 2 in place of (r 1 -) 2 
in Eq.©for B b and use white noise for £f(t): C b {t - t') -> r b S(t - t'). 

The random choice of Xi from trial to trial effectively samples across the 
neuronal populations, so we can then collect the statistics r a , (r°j) 2 (or, equiva- 
lently, (<5r°) 2 ), and C a (t — t') from these trials. These can be used to generate 
an improved estimate of the input noise statistics to be used in 10 in a second 
set of trials, which yields new spike statistics again. This procedure is iterated 
until the input and output statistics agree. This may take up to several hun- 
dred iterations, depending on network parameters and how the computation is 
organized. 

If one tries this procedure in its naive form, i.e., using the output statistics 
directly to generate the input noise at the next step, it will lead to big oscillations 
and not converge. It is necessary to make small corrections (of relative order 
l/-Ko) to the previous input noise statistics to guarantee convergence. 

When one computes statistics from the trials in any iteration, the simplest 
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procedure involves calculating not l(TU|l . but rather Cb(t — t') fEa. lfHfl) '). From 
it, we can proceed in two ways. In the first, from its \t — t'\ — > oo limit we can 
obtain (Sr b ) 2 , and thereby (r b ) 2 = r 2 + (<5r^) 2 for use in calculating Bb in ©. 
Subtracting this limiting value from Cf,(t — t') give us C\(t — t') (which vanishes 
for large \t — t'\) for use in generating the noise £? b (t). This is the first of the 
three methods described above. 

Alternatively, we can use the third method: At each step of our iterative 
procedure we can generate noise directly with the correlations Cb(t — t') (which 
are long-ranged in time) and use r\ in place of (r5) 2 in calculating Bf, (JHJ). 
We have verified that the two methods give the same results when carried out 
numerically, though the second procedure converges more slowly. 

While the true rates in the stationary case are time- independent and C a (t, t') 
is a function only of t — t 1 ', the statistics collected over a finite set of noise- 
driven trials will not exactly have these stationarity properties. Therefore we 
improve the statistics and impose time-translational invariance by averaging 
the measured r a {t) and (<5r"(i)) 2 over t and averaging over the measured values 
C a (t,t') with a fixed t-t'. 

After the iterative procedure converges, so that we have a good estimate 
of the statistics of the input, we want to run many trials on a single neuron 
and compute its firing statistics. This means that the numbers x? b (b = 0, 1, 2J 
should be held constant over these trials. In this case it is necessary to subtract 
out the large t — t' limit of C a (t — t') and use fixed xf b (constant in time and 
across trials) to generate the input noise. (If we did it the other way, without 
the subtraction, we would effectively be assuming that x^ b changed randomly 
from trial to trial, which is not correct.) 

In our calculations we have used 10000 trials to calculate these single-neuron 
firing statistics. We perform the subtraction of the long-time limit of C a (t — t') 
at \t — t'\ = 50, and we have checked that is flat beyond this point in all 
the cases we have done. 

If we perform this kind of measurement separately for many values of the 
if, we will be able to see how the firing statistics vary across the population. 
Here, however, we will confine most of our attention to what we call the "average 
neuron": the one with the average value (0) of all three xf'. 

In particular, we calculate the mean spike count in the 100-ms trials and 
its variance across trials. From this we can get the Fano factor F (the vari- 
ance/mean ratio). We also compute the autocorrelation function, which offers 
a consistency check, since the Fano factor can also be obtained from 

1 Z" 00 

F = - C(r)dr. (18) 

' J —oo 

(This formula is valid when the measurement period is much larger than the 
time over which C(r) falls to zero.) 

We will study below how these firing statistics vary as we change various 
parameters of the model: the input rates ro, parameters that control the balance 
of excitation and inhibition, and the overall strength of the synapses. This will 
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give us some generic understanding of what controls the degree of irregularity 
of the neuronal firing. 



Nonstationary Case 

When the input population is not firing at a constant rate, almost the same 
calculational procedure can be followed, except that one does not average mea- 
sured rates, their fluctuations or correlation function over time. To start out, 
we get initial instantaneous rate estimates from the balance condition, assuming 
that the time-dependent average input currents do not vary too quickly. (This 
condition is not very stringent; van Vreeswijk and Sompolinsky showed that the 
stability eigenvalues are proportional to \/K , so if they have the right sign the 
convergence to the balanced state is very rapid.) 

To do the iterative procedure to satisfy the self-consistency conditions of the 
theory, it is simplest to use the second of the two ways described above (not 
doing any subtraction until the final calculations with single neurons). In this 
case the expression (JHJ does not include rate fluctuations, and we get equations 
for the noise input currents just like (J2J), (@J and © except that the are 
i-dependent and the correlation functions C\, and depend on both t and t' , 
not just their difference. 

The only tricky part is the subtraction of the long-time limit of the correla- 
tion function, which is not simply defined. 

We treat this problem in the following way. We examine the rate-normalized 
quantity 

Da{t,t')= ^ 

r a {t)r a (t') 

We find that this quantity is time-translation invariant (i.e., a function only 
of t — t') to a very good approximation, so we perform the subtraction of the 
long-time limit on it. Then multiplying the subtracted D by r a (t)r a (t') gives 
a good approximation to the true correlation function C a (t,t'). The meaning 
of this finding is, loosely speaking, that when the rates vary (slowly enough) in 
time, the correlation functions just inherit these rates as overall factors without 
changing anything else about the problem. 

We will use the this time-dependent formulation below to simulate exper- 
iments like those of Gershon et al where the LGN input ro(i) to visual 
cortical cells is time-dependent because of the flashing-on and off of the stimu- 
lus. 



3 Results 

The results presented in this chapter were obtained from simulations with pa- 
rameters corresponding to population sizes of N\ = 40,000 excitatory neurons 
and N2 = 10,000 inhibitory neurons. With the above mentioned connection 
probabilities of K a /N a = 0.1, this translates to an average number of K\ = 4000 
excitatory inputs and K2 = 1000 inhibitory inputs to each neuron. The average 
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number of external (excitatory) inputs K was chosen to be equal to K 2 . All 
neurons have the same membrane time constant r of 10 ms. 

To study the effect of various combinations in synaptic strength, we use the 
following generic form to define the intra-cortical weights J a b'- 



For the synaptic strengths from the external population we use J\o — 1 and 
J20 = e. With this notation, g determines the strength of inhibition relative 
to excitation within the network, and e the strength of intracortical excitation. 
Additionally, we scale the overall strength of the synapses with a multiplicative 
scaling factor denoted J s so that each synapse has an actual weight of J s ■ J a b, 
regardless of a and b. 

Figure |21 summarizes how the firing firing statistics depend on all of the 
parameters g, e, and J s . The irregularity of spiking, as measured by the Fano 
factor, depends most sensitively on the overall scaling of the synaptic strength, 
J s . The Fano factor increases systematically as J s increases, and higher values 
of intracortical excitation e also result in higher values of F. The same pattern 
holds for stronger intracortical inhibition, parameterized by g. For all of these 
cases the mean firing rate remains virtually unchanged due to the dynamic 
balance of excitation and inhibition in the network, whereas the fluctuations 
increase with the increase of any of the synaptic weights. 

Interspike interval (ISI) distributions are shown in Figure|3]for three different 
values of J s , keeping e and g fixed at 0.5 and 1, respectively. For a Poisson spike 
train, the Fano factor F = 1, while F > 1 (which we term "superpoissonian" ) 
indicates a tendency of spikes occurring in clusters separated by accordingly 
longer empty intervals, and F < 1 ( "subpoissonian" ) indicates more regularity, 
reflected by a narrower distribution. We have adjusted the input rate ro so that 
the output rate is the same in all three cases. 

The top panel of Figure [3] shows the ISI distribution of a superpoissonian 
spike train, obtained for J s — 1.5. Overlayed on the histogram of ISI counts 
is an exponential curve indicating a Poisson distribution with the same mean 
ISI length. Compared with the Poisson distribution, the superpoissonian spike 
train contains more short intervals, as seen by the peak at short lengths, and 
also more long intervals, causing a long tail. Necessarily, the interval count 
around the average ISI length is lower than that for a Poisson spike train. 

The ISI distribution in the middle panel of Figure [3] belongs to a spike train 
with a Fano factor close to one, obtained for J s = 0.75. The overlayed exponen- 
tial reveals a deviation from the ISI count: while intervals of diminishing length 
are the most likely ones for a real Poisson process, our neuronal spike trains al- 
ways show some refractoriness reflected by a dip at the shortest intervals. (We 
have not used an explicit refractory period in our model. The dip seen here 
simply reflects the fact that it takes a little time for the membrane potential 
distribution to return to its steady-state form after reset.) Apart from this de- 
viation, however, there is a close resemblance between the observed distribution 
and the "predicted" one. 
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Finally, the lower panel of Figure |31 depicts a case with F < 1, with weaker 
synapses, leading to a stronger refractory effect and (since the rate is fixed) an 
accordingly narrower distribution around the average ISI length, as compared 
to the overlayed Poisson distribution. This distribution was obtained with weak 
synapses produced by a small scaling factor of J s = 0.375. 

As mentioned in the previous section, the Fano factor can also be obtained 
by integrating over the spike train autocorrelation divided by the spike rate 
(|18l) . For a Poisson process the autocorrelation vanishes for all lags different 
from zero. In contrast, F > 1 (superpoissonian case) implies a positive integral 
over non-zero lags, whereas in the subpoissonian case there must be a negative 
area under the curve. Figure0]shows examples of autocorrelations for all of the 
three cases. For the superpoissonian case (red dash-dot line), there is a "hill" 
of positive correlations for short intervals, reflecting the tendency toward spike 
clustering. The subpoissonian autocorrelation (blue dotted line) shows a valley 
of negative correlations for short intervals, indicating well separated spikes in a 
more regular spike train. The curve labeled as Poisson (black solid line) does 
have a small valley around zero lag, which reflects once more the refractoriness 
of neurons to fire at extremely short intervals, unlike a completely random 
Poisson process. (Actually, the measured F in this case is slightly greater than 
1, implying that in this case the integral of the very small positive tail for t > 2 
ms is slightly larger than that of the (more obvious) negative short-time dip.) 

Measurements on VI neurons in awake monkeys (see for example Gershon 
et al. QTj) suggest a linear relationship between the log variance and the log 
mean of stimulus-elicited spike counts. We find a similar dependence for neurons 
within our model network. Figure shows results for three different values of 
J s . In each case, five different values of the external input rate ro were tried, 
causing various mean spike counts and variances. The logarithm of the spike 
count variance is plotted as a function of the logarithm of the spike count mean, 
and a solid diagonal line indicates the identity, i.e, a Fano factor of exactly 1. 
We see that for the largest value of J s used here, the data look qualitatively like 
those from experiments, with Fano factors in the range around 1.5 to 2. 

Nonstationary Case 

The results presented in the previous section were obtained with stationary 
inputs, while experimental data like those from |17j were collected from visual 
neurons subject to time-dependent inputs. Therefore, we performed calculations 
of the spike statistics in which the input population rate ro was time-dependent. 
The modeled temporal shape of ro (t) is depicted in Figure El It is the sum of 
three terms: 

r (t) = Ro + A(t) + B(t) (21) 

The first, Ro, is just a constant, as in the preceding section. The second term, 
A(t), rises to a maximum over a 25- ms interval, remains constant for 50 ms, 
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A 


0.375 


0.375 


0.500 


0.500 


0.750 


0.750 


1.000 


1.000 


B 


0.125 


0.375 


0.125 


0.375 


0.250 


0.750 


0.250 


0.750 


F 


1.14 


1.2 


1.22 


1.23 


1.29 


1.36 


1.37 


1.4 


A 


1.500 


1.500 


2.000 


2.000 


3.000 


3.000 


4.000 


4.000 


Bo 


0.500 


1.500 


0.500 


1.500 


1.000 


3.000 


1.000 


3.000 


F 


1.48 


1.5 


1.55 


1.53 


1.57 


1.41 


1.43 


1.34 



Table 1: Stimulus parameters Ao and Bo for the results depicted in Figure 
and the resulting Fano factors F. 

and then falls off to zero over the final 25 ms. 

( 0.5A (1 - cos(4t7r/T)) for < t < T/4 

A(t) = < A for T/4 < t < 3T/4 (22) 

[ 0.5A (1 - cos(4(T - t)-K/T)) for 3T/4 < t < T, 

where T is the total simulation interval of 100 ms. The third term, Bo, rises 
to a maximum in the first 25 ms and then falls back to zero in the next 25 ms, 
remaining zero thereafter. 

( 0.5B (1 - cos(4t7r/T)) for < t < T/4 

B(t) = < 0.5B (1 - cos(4(T/2 - t)ic/T)) for T/4 < t < T/2 (23) 
[ for T/2 < t < T. 

Figure shows the logarithm of the spike count variance plotted against the 
logarithm of the spike count mean for various non-stationary inputs character- 
ized by different values of A and B - The graph shows results for J s = 1, 
e = 0.5, g = 1, and a background rate of R = 0.1. Table^shows the choice of 
the sixteen combinations of the stimulus parameters Ao and Bo, together with 
the resulting Fano factors F for the simulated neuron. 

The data look qualitatively like those obtained from in-vivo experiments 
|17| and are similar to the superpoissonian case in Figure The neuron fires 
consistently in a superpoissonian regime with Fano factors slightly higher than 
1 and an almost linear relationship between the log variance and the log mean 
for low spike counts. For higher spike counts, the curve bends towards values of 
lower Fano factors, just as for stationary inputs (Figure |SJ. In both cases, this 
bend reflects the the decrease in irregularity of firing caused by an increasingly 
prominent role of refractoriness for shorter interspike intervals. 

4 Discussion 

Cortical neurons receive thousands of both excitatory and inhibitory inputs, 
and despite the high number of inputs from nearby neurons with similar firing 
statistics and similar connectivity, their observed firing is very irregular [5J 1101 
[IIl[T2[I3[Il[T3[Tl[n|[Tl[T3|iniin]- Dynamically balanced excitation and 
inhibition through a simple feedback mechanism provides an explanation that 
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naturally accounts for this phenomenon without requiring fine tuning of the 
parameters |2 03 GI 03 • Moreover, neurons in such model networks show an 
almost linear input-output relationship (input current versus firing frequency), 
as do neurons in the neocortex. 

Here, we have extended the mean-field description of the dynamically bal- 
anced asynchronous firing state to analyze firing correlations. We found that 
the relationship between the observed irregularity of firing (spike count vari- 
ance) and the firing rate (spike count mean) of the neurons resemble closely 
data collected from in-vivo experiments (see Figures and . To do this, we 
developed a complete mean-field theory for a network of leaky integrate-and- 
fire neurons, in which both firing rates and correlation functions are determined 
self-consistently. Using an algorithm that allows us to find the solutions to 
the mean-field equations numerically, we could elucidate how the strength of 
synapses within the network influences the expected firing statistics of cortical 
neurons in a systematic manner (see Figure EJ. 

We have shown that the irregularity of firing, as measured by the Fano 
factor, increases with increasing synaptic strengths (Figure |2J). Nearly Poisson 
statistics (with F w 1) are observed for moderately strong strengths, but the 
transition from subpoissonian to superpoissonian statistics is smooth, without 
a special role for F = 1. 

The higher irregularity in the spike counts is always accompanied by a ten- 
dency toward more "bursty" firing. (These bursts are a network effect; the 
model contains only leaky integrate-and-fire neurons, which do not burst on 
their own.) This burstiness can best be seen in the spike train autocorrelation 
function (Figure EJ, which acquires a hill of growing size and width around zero 
lag for increasing Fano factors. The interdependence between firing irregular- 
ity and bursting can be understood with help of the ISI distributions depicted 
in Figure |21 when the rate, and thus the average ISI, is kept constant, then 
any higher count for shorter-than-average ISIs must be accompanied by an ac- 
cordingly higher count for longer ISIs (indicating bursts), and vice versa. Thus 
higher irregularity always goes hand in hand with a higher tendency toward 
temporal clustering of spikes. 

Why do stronger synapses lead to higher irregularity in firing? The size of 
the input current fluctuations in JBJ are controlled by the J a b, and so, therefore, 
are the corresponding membrane potential fluctuations. Thus, for example, the 
width of the steady-state membrane potential distribution is proportional to 
J s . We next have to consider where this distribution is centered. Remembering 
that, according to the balance condition, the firing rate is independent of J s , 
the center of the distribution has to move farther away from threshold as J s is 
increased in order to keep the rate fixed. Therefore, for very small J s almost 
the entire equilibrium membrane potential distribution will lie well above the 
post-spike reset value, while for large J s it will be mostly below reset. 

Immediately after a spike, the membrane potential distribution is a delta- 
function centered at the reset (here 0). It then spreads and its mean moves 
up or down toward its equilibrium value. This equilibration will take about 
a membrane time constant. If the equilibrium value is well above zero (the 
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small- J s case), the probability of reaching threshold will be suppressed during 
this time, implying a refractory dip in the ISI distribution and the correlation 
function and a tendency toward a Fano factor less than 1. 

In the large- J s case, on the other hand, where the membrane potential is 
reset much closer to the threshold than to its eventual equilibrium value, the 
initial rapid spread (with the width growing proportional to J s Vi) leads to 
an enhanced probability of early spikes. At short times this diffusive spread 
dominates the downward drift of the mean (which is only linear in t). Thus 
there is extra weight in the ISI distribution and a positive correlation function 
at these short times, leading to a Fano factor greater than 1. 

Empirically, an approximate power-law relationship between the mean and 
variance of the spike count has frequently been observed for cortical neurons (see, 
e -g-> [TT1IT5HT7II5D] '). Our model shows the same qualitative feature (Figures 
andEJ), though we have no argument that the relation should be an exact power 
law. However, this agreement suggests that the model captures at least part of 
physics underlying the firing statistics. 

As already observed, not all of the variability in measured neuron responses 
has to be explained in the manner outlined above. Changing conditions dur- 
ing the run of a single experiment may introduce extra irregularity, caused by 
collecting statistics over trials with different mean firing rates. The present anal- 
ysis shows why - and how much - irregularity can be expected due to intrinsic 
cortical dynamics. 

Our formulation of the mean-field theory is general enough to allow straight- 
forward extensions to greater biological realism and to more complicated net- 
work architectures. We have introduced a generalization of this model with 
conductance-based synapses in another paper [21]. We have also extended the 
model to include systematic structure in the connections, modeling an orienta- 
tion hypercolumn in the primary visual cortex |25j . Moreover, our algorithm for 
finding the mean-field solutions is not restricted to networks of integrate-and- 
fire neurons. It can be applied to any kind of neuronal model. Furthermore, 
any kind of synaptic dynamics can be incorporated by using synaptically filtered 
spike trains to compute the self-consistent solutions. 
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Figure 2: Fano factors as a function of overall synaptic strength J s 
and intracortical excitation strength e for three different inhibition 
factors: g = 1, 1.5, and 2, respectively. The increase of any 
of these parameters results in more irregular firing statistics as 
measured by the Fano factor. 
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Figure 3: Interspike interval distributions for fixed e = 0.5 
and g = 1, and three different values of overall synaptic 
strength J s : 1.5 (superpoissonian), 0.75 (Poissonian), and 
0.375 (subpoissonian). Overlayed on each figure is the expo- 
nential fall-off of a true Poisson distribution with the same 
average rate as in all of the three cases. 
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Figure 4: Three different spike train autocorrelations illustrating the relation- 
ship between the Fano factor F and the area under the curve. For F = 1 (Pois- 
sonian, black solid line), the autocorrelation is an almost perfect delta function. 
F > 1 (superpoissonian, red dash-dot line) is reflected by a hill generating a 
positive area, and F < 1 (subpoissonian, blue dotted line) is accompanied by a 
valley of negative correlations. (See the text for more details.) 
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Figure 5: Spike count log(variance) vs. log(mean) for three different values 
of overall synaptic strength J s , varying the external input rate ro- For J s = 
1.25 (superpoissonian, red triangles) the data look qualitatively like those from 
experiments. The other values for J s are 0.75 (Poisson, black stars) and 0.375 
(subpoissonian, blue crosses). 
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Figure 6: Parametrization of the time-dependent input rate ro(t). The input is 
modeled as the sum of three functions: (1) a stationary background rate (which 
is zero in this case); (2) a tonic part, which rises within the first 20 ms to a 
constant level of Aq where it stays for 60 ms, falling back to zero within the last 
20 ms; and (3) an initial phasic part, which is nonzero only in the first 50 ms, 
rising to a maximum value of Bq. 
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Figure 7: Spike count log(variance) vs. log(mean) for time-varying external 
inputs with varying overall strength. The neuron in the simulated network (red 
triangles) fires in a superpoissonian regime, with an almost linear relationship 
for low spike rates between the log variance and the log mean, resembling closely 
data obtained from in- vivo experiments. The diagonal solid line indicates the 
identity of variance and mean (Fano factor F = 1). 
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